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introduction: 


Under  a  prior  research  effort  entitled  "Flow  around  Maneuvering  Appended  Bodies'"  the 
NSF  Engineering  Research  Center  at  Mississippi  State  University  (MSU/ERC)  and  the 
Applied  Research  Laboratory  at  The  Pennsylvania  State  University  (ARL)  developed  a 
physics-based  means  of  predicting  the  maneuvering  characteristics  of  a  self-propelled, 
fully-appended  underwater  vehicles.  Subsequently,  a  follow-on  project  to  extend  the  initial 
capability  was  undertaken.  This  effort  addresses  many  of  the  issues  that  remained  at  the 
conclusion  of  the  original  project.  Specific  areas  being  addressed  include:  (1) 
computation  of  flow  at  very  high  (full-scale)  Reynolds  number,  (2)  investigation  of  different 
grid  topologies,  (3)  turbulence  models  studies,  (4)  effective  use  of  parallel  processing,  and 
(5)  Depth-changing  maneuvers.  Subsequently,  the  follow-on  effort  was  divided  into  two 
parts  to  allow  consideration  of  classified  configurations  under  a  separate  project. 


1 .  Computation  of  flow  at  very  high  (full-scale)  Reynolds  number: 

To  evaluate  the  issues  involved  in  conducting  high  Reynolds  number  computations,  a 
case  study  is  set  up  to  perform  Navier-Stokes  computations  at  a  very  high  Reynolds 
number. 

SUBOFF,  a  model  scale  submarine,  has  been  one  of  the  geometries  used  for  the 
validation  of  the  method  developed  in  the  course  of  the  above  mentioned  research  effort. 
The  Reynolds  number  is  12  million  for  the  experimental  and  computational  studies 
performed  with  SUBOFF  (based  on  body  length).  However,  the  full  scale  Reynolds 
number  for  submarines  is  about  two  orders  of  magnitudes  higher  than  the  one  used  for 
the  model  scale.  Therefore,  in  the  case  study,  a  Reynolds  number  of  1.2  billion  (based  on 
body  length)  was  used,  corresponding  to  a  full-scale  submarine  operating  Reynolds 
number.  The  SUBOFF  geometry  with  four  stern  appendages  used  for  this  case  is  shown 
in  Figure  1 .  Since  a  grid  for  this  case  already  existed,  it  was  re-gridded  to  satisfy  the  high 
Reynolds  requirement  for  resolving  the  boundary  layer. 

The  near-body,  high  grid  resolution  requirement  may  increase  the  degree  of  difficulty  for 
the  timely  creation  of  the  grids  with  the  currently  available  grid  generation  packages.  The 
cell  aspect  ratio  increases  substantially  and  generally  it  is  much  harder  to  create  a  good 
quality  grid.  This  tests  the  robustness  of  the  flow  solver  and  may  increase  the  number  of 
cycles  required  for  convergence.  Accurate  definition  of  the  underlying  body  geometry 
(with  CAD)  is  also  a  difficult  requirement.  Typically  the  manufacturing  tolerance 
requirements  for  the  geometry  are  much  larger  than  the  high  Reynolds  number  grid 
tolerance  requirements,  and  many  CAD  programs  are  not  well-suited  to  providing  the  high 
tolerance  needed  for  RANS  grid  generation. 

Figure  2  shows  the  comparison  between  the  computed  velocity  profile  (using  the  Baldwin- 
Lomax  turbulence  model)  and  the  law  of  the  wall  for  the  Reynolds  number  of  1 2  million 
(based  on  the  body  length).  This  velocity  profile  is  taken  on  the  parallel  mid-section  of  the 
hull.  It  can  be  seen  that  the  velocity  profile  has  been  predicted  with  approximately  two 


grid  points  in  the  viscous  sublayer.  Figure  3  shows  the  velocity  profile  at  the  same  axial 
location  for  the  Reynolds  number  computations  of  1 .2  billion.  The  total  number  of  grid 
points  used  for  both  computations,  as  well  as  grid  density  on  the  hull  and  appendages, 
are  also  shown  in  these  figures.  The  velocity  profile  for  the  high  Reynolds  number 
computation,  normal  to  an  appendage  at  r/R  =  0.8,  where  r  is  the  radial  distance  from  the 
centerline  and  R  is  the  radius  of  the  hull,  is  shown  in  Figure  4. 

Figures  5  and  6  show  the  integrated  axial  force  coefficient  for  both  computations.  The 
contribution  of  the  viscous  stresses  and  pressure  to  the  axial  force  are  also  shown  in 
these  figures.  It  is  noted  that  the  number  of  cycles  required  to  convergence  has  increased 
dramatically  for  the  higher  Reynolds  number  case.  A  better  quality  grid  can  decrease  the 
number  of  cycles  required  to  achieve  convergence.  This  was  demonstrated  by  MSU/ERC 
in  a  similar  computation  with  an  improved  grid.  The  need  to  create  good  quality  grids  in  a 
timely  manner  lead  to  the  following  study. 


2.  Investigation  of  different  grid  topologies 

A  study  to  improve  grid  quality  and  simplify  grid  generation  for  complex  geometries  was 
initiated  under  a  separate  propulsor  analysis  effort.  The  lessons  learned  from  that  study 
were  applied  to  a  submarine  configuration  as  described  below. 

This  study  intended  to  investigate  the  implementation  of  alternate  grid  topologies,  to 
explore  ways  to  accelerate  the  grid  generation  process  and  to  improve  the  quality  of  the 
final  grid.  The  geometry  used  for  this  study  is  SUBOFF  with  four  stern  appendages.  Most 
early  versions  of  the  flow  solver  UNCLE  require  a  certain  type  of  grid  topology  namely, 
CH-grid  type.  This  requirement  imposes  restrictions  on  the  quality  of  the  grid  and  on  the 
density  of  the  grid  points  near  solid  surfaces;  on  certain  geometries,  an  0-grid  or  an  OH- 
grid  may  yield  a  better  quality  grid.  More  recent  versions  of  the  UNCLE  code  do  not  have 
this  grid  topology  restriction. 

For  this  study,  two  grid  generation  software  packages  were  considered:  EAGLEVIEW  and 
GRIDPRO.  EAGLEVIEW  was  used  to  create  an  H-grid  around  the  stem  appendages 
whereas  GRIDPRO  was  used  to  create  OH-grid  around  the  hut!  and  appendages.  The 
EAGLEVIEW  grid  was  generated  several  years  ago  by  MSU/ERC.  The  commercially 
available  GRIDGEN  package  is  now  used  by  ARL,  since  EAGLEVIEW  is  no  longer 
supported. 

GRIDPRO  is  designed  to  create  0-grids  around  surfaces.  The  0-grid  may  then  be 
surrounded  by  an  H-grid  configuration.  In  GRIDPRO,  once  the  grid  topology  is  defined, 
the  grid  generation  is  automatic.  However,  if  the  topology  is  not  satisfactory  GRIDPRO 
will  not  converge  and  therefore,  will  not  yield  a  usable  grid.  GRIDPRO  creates  a  large 
number  of  blocks,  which  with  the  help  of  a  postprocessor  can  be  reduced.  These  blocks 
are  located  somewhat  arbitrarily  relative  to  one  another  with  arbitrary  alignment  of  the 
computational  coordinate  directions  that  makes  it  difficult  to  distinguish  the  block 
interfaces.  Therefore,  a  preprocessor  was  developed  to  identify  the  inter-block  surfaces 


and  set  up  an  important  input  file  for  the  parallei  code.  Figure  7  shows  the  overall 
blocking  structure  used  by  EAGLEVIEW  to  create  a  CH-grid  and  the  blocking  structure 
created  by  GRIDPRO  to  generate  the  OH-grid.  It  should  be  noted  that  the  block  structure 
created  by  GRIDPRO  is  automatically  created  by  the  software  from  the  user-defined 
topology,  whereas  with  EAGELVIEW  the  user  explicitly  defines  the  structure.  Figures  9 
and  10  show  the  grid  around  an  appendage  in  the  stern  region  for  both  the  H-grid  and  O- 
grid.  It  should  be  noted  that  the  grid  points  around  the  appendage  are  much  denser  for 
the  OH-grid  than  the  H-grid  (by  approximately  one  order  of  magnitude).  The  first  grid 
point  off  the  appendage  surface  is  1.1  x  10'®  for  the  OH-grid  and  4.9  x  10‘®forthe  H-grid, 
both  in  units  of  body  length. 

Three  computations  were  performed  to  examine  the  grid  created  by  the  GRIDPRO  for 
SUBOFF  with  four  stern  appendages.  These  computations  are  performed  at  0,  8  and  18 
degree  angle  of  drift.  Velocity  profiles,  taken  at  the  mid-section  of  the  hull,  for  both  grid 
types  are  shown  in  Figures  2  and  11.  It  should  be  noted  that  the  distribution  of  the  grid 
points  normal  to  the  hull  is  much  denser  for  the  OH-grid  than  the  H-grid.  The  viscous 
sublayer  is  defined  with  eight  grid  points  for  the  OH-grid  whereas  it  is  defined  with  only 
two  grid  points  for  the  H-grid.  This  high  resolution  with  the  OH-grid  guarantees  that  the 
viscous  stresses  are  well  predicted. 

The  convergence  history  for  the  axial  force  coefficient  is  shown  in  Figure  12.  Both  cases 
were  run  with  the  Baldwin-Lomax  model.  As  shown,  the  OH-grid  takes  longer  to  converge 
than  H-grid.  This  is  partly  because  the  multi-grid  scheme  could  not  be  used  for  the  OH- 
grid,  since  the  number  of  grid  points  in  each  block  was  not  set  to  allow  multi-grid.  Also, 
the  higher  cel!  aspect  ratio  and  near-surface  resolution  might  reduce  the  convergence 
rate,  especially  in  the  absence  of  multi-grid  iterations. 

Further,  the  Baldwin-Lomax  model  with  a  search  parameter  of  2500  was  used  for  the  flow 
computation  for  the  H-grid.  The  search  parameter  use  is  described  in  the  attached 
publication.  Appendix  A.  No  search  parameter  was  specified  for  the  OH-grid  computation, 
so  that  the  search  for  Fmax  extended  to  the  far  field  boundary. 

The  variation  of  the  axial  force  coefficient  versus  angle  of  drift  for  the  two  grid  topologies, 
CH  and  OH,  are  shown  in  Figure  13.  The  experimental  data  is  also  shown  in  this  figure. 
For  small  angles  of  drift,  the  computation  performed  on  the  OH-grid  is  closer  to  the 
experimental  data  than  the  CH-grid  computations.  For  angles  of  drift  between  10  and  18 
degrees,  the  CH-grid  results  are  closer  to  the  experimental  results  for  negative  values  of 
angle  of  drift.  It  should  be  noted  that  the  original  Baldwin-Lomax  turbulence  model  is  used 
for  the  latter  two  calculations.  A  prediction  using  Coakley’s  q-co  turbulence  model  is  also 
shown  in  the  figure.  Figure  13  also  shows  the  prediction  for  CH-grid  using  Baldwin-Lomax 
turbulence  model  with  the  search  parameter  of  2500.  A  very  good  comparison  with  the 
experimental  data  is  observed  for  positive  values  of  the  angle  of  drift.  Because  of  the 
asymmetry  in  the  experimental  data  (for  a  symmetric  body),  no  model  can  predict  these 
experimental  data  for  both  positive  and  negative  values  of  the  drift  angle.  The  asymmetry 
in  the  data  indicates  an  additional  error  not  captured  by  the  uncertainty  bars  shown.  A 
more  detailed  analysis  of  the  choice  of  the  turbulence  models  on  the  predicted  results  is 


given  in  the  enclosed  publication.  The  present  results  indicate  the  need  for  a  systematic 
grid  resolution  study  to  properly  quantify  the  grid-related  errors.  This  must  be  done  before 
a  more  thorough  analysis  of  the  turbulence  model  can  be  performed. 

Comparison  of  the  predicted  pressure  coefficients  and  the  experimental  data  are  shown  in 
Figures  14  and  15.  The  only  deviation  between  the  experimental  data  and  the 
computational  results  occurs  at  the  tip  of  the  appendage  for  the  H-grid.  This  is 
understandable,  since  the  tip  of  the  appendage  in  this  grid  is  "pinched"  and  the  tip  is 
shaped  like  a  triangular  prism,  rather  than  being  flat,  which  is  the  case  with  the  OH-grid. 
Apparently,  this  did  not  lead  to  significant  errors  in  the  overall  force  and  moment 
predictions  for  the  captive  model  tests.  Note  that  cap  blocks  are  required  for  either 
topology  to  accurately  capture  the  actual  (flat)  tip  geometry. 

Comparison  of  the  predicted  skin  friction  for  the  OH-grid  with  the  experimental  data  is 
shown  in  Figure  16.  Considering  the  magnitude  of  the  skin  friction,  the  prediction  is  very 
reasonable. 

An  important  feature  of  a  grid  is  its  ability  to  capture  flow  features  of  interest  in  the 
computational  domain.  Figure  9  shows  grid  distribution  around  an  appendage,  using  the 
H-grid  created  with  EAGLEVIEW  and  the  OH-grid  created  with  GRIDPRO.  As  shown  in 
this  figure,  grid  points  for  the  H-grid  are  more  densely  distributed  downstream  of  the 
trailing  edge.  This  OH-grid  was  not  constructed  to  capture  the  appendage  wake; 
therefore  the  H-grid  solution  better  captures  the  wake  and  the  flow  aft  of  the  trailing  edge. 
Figures  17  and  18  show  the  axial-component  of  the  velocity  on  the  appendages  at  18- 
degree  angle  of  drift  for  the  H-grid  and  for  the  OH-grid,  respectively.  The  H-grid 
qualitatively  predicts  expected  flow  features  such  as  shed  vortices.  Considering  cases 
where  a  propulsor  exists  immediately  downstream  of  the  fins,  these  vortices  will  interact 
with  the  propulsor  and  will  change  the  loading  on  the  propulsor.  As  indicated  by  these 
figures  this  particular  OH-grid  fails  to  capture  the  flow  features  shown  by  the  H-grid  aft  of 
the  trailing  edge.  However,  in  the  absence  of  experimental  flow  data  in  this  region,  it  is 
not  possible  to  make  a  definitive  statement  about  the  required  resolution.  This 
demonstrates  the  need  for  further  work  to  quantify  the  required  resolution  of  salient  flow 
features  needed  to  accurately  predict  propulsor  side  forces  during  maneuvering  situations. 

Figures  19(a)-(b)  show  the  pressure  distribution  on  the  SUBOFF  for  both  grid  topologies 
at  18-degree  incidence  angle.  Contours  of  pressure  on  a  plane  perpendicular  to  the  body 
are  also  shown  in  these  figures.  It  should  be  noted  that  in  the  H-grid  case,  the  Baldwin- 
Lomax  turbulence  model  with  the  search  parameter  of  2500  was  used  and  in  the  case  of 
the  OH-grid,  the  search  parameter  is  not  limited.  To  capture  the  vortical  structure  as 
depicted  by  the  H-grid  computation,  it  might  be  necessary  to  refine  the  OH-grid  or  to 
implement  the  search  parameter  limitation.  The  computed  pressure  distribution  in  the 
stern  region  at  the  18-degree  drift  angle  for  the  OH-grid  is  shown  in  Figure  20. 


3.  Turbulence  model  investigation: 

A  modified  version  of  the  Baldwin-Lomax  model  has  been  developed,  and  various 
versions  of  the  k-e  and  q-co  models  were  evaluated. 

A  series  of  computations  were  performed  to  evaluate  the  suitability  of  several  algebraic  as 
well  as  two-equation  turbulence  models.  The  algebraic  models  considered  are  the 
original  Baldwin-Lomax  model,  the  Degani  and  Schiff  variation  of  the  Baldwin-Lomax 
model,  and  a  new  proposed  model.  Furthermore,  the  two-equation  model  of  Coakley  was 
also  evaluated.  For  the  CH-grid  used,  the  computations  demonstrate  that  the  original 
Baldwin-Lomax  is  not  suitable  for  the  three-dimensional  computations  that  include  vortical 
structure  away  from  the  body.  This  model  dissipates  and  distorts  the  vortical  structures 
and  underestimates  axial  and  lateral  forces.  The  modified  Baldwin-Lomax  model  with  the 
specified  search  parameter  produces  well-defined  vortices,  and  the  computed  forces 
compare  well  with  the  experimental  data  if  a  suitable  search  parameter  is  used.  However, 
obtaining  the  best  value  of  the  search  parameter  relies  on  trial  and  error;  therefore,  these 
models  are  not  general. 

A  new  proposed  algebraic  model  is  a  more  general  mode!  that  does  not  need  a  search 
parameter  and  does  not  require  prior  knowledge  of  the  location  of  the  leeward  or 
windward  side  of  the  hull.  It  preserves  the  structure  of  the  vortices  better  than  the  original 
Baldwin-Lomax  model,  although  not  as  good  as  when  the  optimal  search  parameter  is 
specified.  However,  in  its  present  form,  it  assumes  that  the  core  of  the  vortical  structures 
is  aligned  with  the  centerline  axis  of  the  body.  A  more  general  extension  of  the  model 
should  consider  alignment  of  the  vortical  structures  in  any  direction.  The  two-equation 
model  of  Coakley  was  also  considered  for  several  computations.  This  model  avoids  any 
problem-specific  tailoring  and  is  able  to  predict  the  forces  on  the  body;  however,  a  price  is 
paid  in  terms  of  efficiency  and  robustness.  In  addition,  the  eddy  viscosity  appears  to  be 
too  high  in  the  cores  of  the  shed  vortices,  and  so  the  vortices  diffuse  too  quickly. 

More  than  thirty  cases  were  computed  to  study  various  versions  of  the  Baldwin-Lomax 
model  and  two-equation  models.  A  detailed  description  of  the  cases  computed  and  the 
results  are  given  in  the  attached  paper.  Appendix  A. 


4.  Parallel  Processing: 

A  major  issue  in  developing  parallel  codes  relates  to  load  balancing  among  the  various 
processors.  The  parallel  UNCLE  code  is  structured  to  run  one  block  per  processor. 
Therefore,  proper  load  balancing  occurs  only  if  all  the  blocks  in  the  domain  have  the  same 
number  of  grid  points.  Under  this  constraint  each  processor  experiences  the  same 
execution  load.  In  practice,  this  assumption  has  some  shortcomings.  First,  it  is  not 
always  practical  or  even  possible  to  split  the  computational  domain  into  equally  sized 
blocks.  This  is  especially  true  when  the  geometry  is  complex  and  contains  numerous 
appendages,  in  these  situations  the  geometry  dictates  breaking  the  computational 
domain  into  blocks  of  varying  sizes,  resulting  in  a  very  unbalanced  loading  of  the 


processors.  Second,  v^hen  the  number  of  blocks  exceeds  the  number  of  available 
processors,  it  will  become  necessary  to  submit  more  than  one  block  to  some  or  all 
processors.  In  order  to  overcome  these  shortcomirngs,  a  'relative  frame’  version  of  the 
parallel  version  of  UNCLE  (denoted  UNCLE-REL)  was  modified  under  another  project  to 
enable  submission  of  any  number  of  blocks  to  a  single  processor.  In  this  way,  balanced 
loading  can  be  obtained  by  bundling  a  collection  of  small  sized  blocks  together  and 
submitting  them  to  a  single  processor.  The  total  number  of  grid  points  in  those  blocks  is 
constrained  only  by  the  available  memory  for  each  processor.  Furthermore,  in  situations 
where  the  number  of  processors  are  less  than  the  number  of  the  blocks,  the  computation 
still  can  be  executed,  provided  that  sufficient  memory  is  available. 

For  verification  purposes  a  calculation  was  performed  with  the  modified  parallel  UNCLE- 
REL  on  a  flat  plate  with  four  blocks  -  each  set  of  two  blocks  are  submitted  to  a  single 
processor.  The  results  are  identical  with  the  results  obtained  using  an  ARL  serial  version 
of  the  UNCLE  code.  In  the  future,  it  would  be  desirable  to  apply  the  same  modifications  to^ 
the  more  general  'clicking'  version  of  the  UNCLE  code,  recently  obtained  from  MSU/ERC. 
This  version  has  been  rewritten  using  many  features  of  FORTRAN  90  and  it  is 
substantially  different  from  UNCLE-REL. 


5.  Depth-changing  Maneuvers: 

Two  maneuvering  computations  are  performed  to  investigate  the  depth-changing 
maneuvers  of  a  fully  appended  submarine  type  vehicle.  The  geometry  used  for  these 
computations  is  that  of  the  SUBOFF.  It  includes:  sail,  sailplanes,  rudders,  and  horizontal 
stern  planes  and  an  open,  five  bladed  propulsor.  It  should  be  noted  that  the  original 
SUBOFF  does  not  include  sailplanes,  nor  does  it  include  a  propulsor.  These  components 
are  added  to  the  original  SUBOFF  to  (a)  make  the  vehicle  a  self-propelled  vehicle  and  (b) 
to  resemble  a  real  submarine  with  all  the  appendages.  However,  there  are  no 
experimental  data  available  to  compare  with  the  computational  simulations,  but  the 
purpose  behind  these  calculations  are  to  demonstrate  the  applicability  of  the  method  and 
to  observe  that  the  characteristics  of  the  maneuvering  simulations  correlates  to  the 
physics  of  the  problem. 

To  induce  a  change  on  the  elevation  of  a  submarine,  typically  the  horizontal  stern  planes 
are  deflected.  To  simulate  the  effects  of  the  deflections  of  the  sternplanes,  an  external 
vertical  thrust  is  applied  on  the  planes.  Based  on  the  magnitude  and  direction  of  these 
forces,  the  integrated  value  of  the  forces  and  moments  on  the  whole  vehicle  changes. 

The  net  result  leads  to  change  in  the  direction  of  motion  of  the  vehicle  and  eventually  to 
the  change  in  the  elevation  of  the  vehicle. 

All  maneuvering  computations  start  from  a  periodic  solution  at  which  only  the  governing 
equations  of  the  fluid  flow  are  numerically  solved.  Once  the  periodic  solution  is  achieved 
equations  governing  vehicle  dynamics  are  iteratively  solved  with  the  flow  equations  to 
simulate  the  maneuver. 


Two  depth-changing  maneuvers  are  performed,  a  rising  maneuvers  and  a  diving 
maneuvers. 

Rising  Maneuver: 

The  trajectory  of  the  vehicle  during  the  rising  maneuver  is  shown  in  Fig.  21.  During  this 
maneuver,  the  vehicle  goes  through  five  distinct  phases.  (1)  Straight-and-level  phase  that 
extends  from  t  =  0.0,  to  t  =  3.0.  During  this  period,  no  external  force  is  applied  to  the 
vehicle  and  it  is  moving  forward  under  the  thrust  induced  by  its  own  propulsor.  (2)  Pitch- 
up  phase  extends  from  approximately  t  =  3.0  to  t  =  4.0.  To  induce  the  pitch-up  motion  on 
the  vehicle,  an  external  force  is  applied  on  the  horizontal  stern  planes.  The  magnitude  of 
the  force  and  its  duration  is  shown  in  the  Fig.  22.  (3)  Straight-and-ievei  flight  phase  that 
extends  from  approximately  t  =  4.0  to  t  =  4.8.  The  force  applied  on  the  stern  plane  in  the 
phase  2  is  removed,  allowing  the  vehicle  to  return  to  the  straight-and-level  flight.  It  should 
be  noted  that  since  the  orientation  of  the  vehicle  has  been  changed  during  the  previous 
pitch-up  phase,  the  vehicle  will  continue  to  move  to  higher  elevations.  (4)  Pitch-down 
phase  starts  from  approximately  t  =  4.8  to  t  =  5.2.  At  this  phase  of  the  maneuver  the 
thrust  on  the  horizontal  stern  planes  are  changed  to  induce  a  pitch-down  maneuver. 

Pressure  distribution  on  the  body  of  the  vehicle  is  shown  in  the  Figs.  21  and  23.  It  is 
worth  noting  the  position  of  the  stagnation  point  on  the  nose  of  the  vehicle  and  compare  it 
with  its  position  in  the  pitch-up  phase  of  the  maneuver,  in  the  pitch-up  phase  the  location 
of  the  stagnation  point  is  on  the  upper  part  of  the  nose,  as  expected,  in  the  pitch-down 
phase  of  the  maneuver,  the  location  of  the  stagnation  point  moves  to  the  opposite  side  on 
the  lower  part,  again  as  physics  of  the  problem  dictates.  (5)  In  the  final  phase  the  force 
applied  on  the  stern  plane  in  the  phase  4  is  removed,  allowing  the  vehicle  to  return  to  the 
straight-and-level  flight.  From  t  =  5.2  to  the  end  of  the  computation  the  vehicle  travels  in 
the  straight-and-level  mode  again. 

Fig.  22  shows  the  magnitude  and  the  manner  in  which  the  appiied  force  on  the  horizontal 
stern  planes  are  changed.  The  trajectory  of  the  location  of  the  vehicle  during  the 
maneuver  is  also  shown  in  this  figure.  Angular  rates  experienced  by  the  vehicle  during 
this  maneuver  are  shown  in  Fig.  24.  The  angles  of  roll,  pitch  and  yaw  are  shown  in  the 
Fig.  25. 

Diving  Maneuver; 

The  procedure  used  to  simulate  the  diving  maneuver  is  very  similar  to  the  procedure  used 
during  the  rising  maneuver,  except  the  direction  of  the  applied  force  on  the  horizontal 
stern  plane. 

The  trajectory  of  the  location  of  the  vehicle  and  the  value  of  the  force  applied  on  the  sten 
plane  is  given  in  the  Fig.  26.  The  angles  and  the  angular  rate  experienced  by  the  vehicle 
are  given  in  the  Figs.  27  and  28.  Considering  the  behavior  of  the  roll  rate  in  both  rising 
and  diving  maneuvers,  it  is  interesting  to  note  that  the  vehicle  experiences  a  far  more 
radical  roll  from  side  to  side  in  the  rising  maneuver  than  in  the  diving  maneuver.  This 


could  be  attributed  to  the  way  that  the  sailplanes  interact  with  the  flow  in  these 
maneuvers.  In  the  pitch-up  maneuver,  sailplanes  directly  interact  with  the  flow  whereas  in 
the  diving  maneuver,  the  hull  obstructs  the  free  flow  of  the  liquid  towards  the  sailplane  and 
thereby  inducing  a  far  lesser  fluctuation  in  the  roll  rate. 

Conclusions: 

A  number  of  issues  pertaining  to  various  aspects  of  the  maneuvering  simulations  are 
addressed.  These  issues  involved  the  choice  for  topology,  the  turbulence  model,  coupling 
of  the  Navier-stokes  equations  of  flow  and  the  six-degree-of-freedom  equations  of  motion, 
full  scale  capabilities  and  the  applicability  of  the  method  to  simulate  real  world 
maneuvering  scenarios. 

The  investigation  of  grid  topologies  demonstrated  that  great  care  must  be  taken  in 
constructing  a  grid  to  preserve  wakes  and  vorticies  downstream  of  appendages.  The  H- 
grid  topology  allows  concentration  of  grid  cells  across  an  appendage  wake  in  a  relatively 
straightforward  manner  when  the  wake  is  closely  aligned  with  the  streamwise  direction. 
The  OH-grid  topology  requires  more  careful  control  of  the  0-bIock  around  an  appendage 
and  the  surrounding  H-block  to  achieve  satisfactory  wake  resolution.  The  0-block  around 
an  appendage  allows  higher  wall-normal  resolution  to  be  achieved  with  better  grid  quality 
(cell  orthogonality)  than  with  an  H-block.  The  proper  prediction  of  the  interaction  of  wakes 
and  vorticies  with  a  propulsor  is  especially  important  in  determining  the  steady  and 
unsteady  propulsor  loads,  which  are  very  important  in  correctly  predicting  the 
maneuvering  behavior  of  a  vehicle. 

Several  turbulence  models  were  investigated.  They  included  the  Baldwin-Lomax 
algebraic  model  and  elliptic  two-equation  models.  The  study  demonstrated  that  the 
original  Baldwin-Lomax  model  is  not  suitable  for  the  three-dimensional  computations  that 
include  vortical  structures  away  from  the  body.  This  model  dissipates  and  distorts  the 
vortical  structures  and  underestimates  axial  and  lateral  forces.  The  algebraic  model  was 
modified  to  account  for  the  effects  of  the  leeward-side  vortices  on  the  underlying  viscous 
layers.  The  modified  Balwin-Lomax  with  an  appropriately  specified  search  parameter 
produces  very  well-defined  vortices,  and  the  computed  forces  compare  well  with  the 
experimental  data  only  if  a  suitable  search  parameter  is  used.  Furthermore,  the  modified 
Baldwin-Lomax  preserved  the  vorticity  better  than  any  of  the  two-equation  models  studied. 
This  conclusion  is  in  agreement  with  the  results  obtained  by  other  investigators. 

Computations  are  presented  to  demonstrate  the  capability  of  the  flow  solver  to  simulate 
very  high  Reynolds  number  (i.e.,  1.2  x  10^)  flow  over  a  submarine-like  configuration.  The 
computation  demonstrated  that  the  very  high  Reynolds  number  computation  can  be 
performed;  however,  the  creation  of  the  high-resolution  grid  in  this  case  presented  a 
challenge. 

Two  depth-changing  maneuvering  simulations  were  performed  to  demonstrate  the 
suitability  of  the  method.  These  simulations,  although  not  compared  with  experiments, 
demonstrated  the  physical  behavior  expected  from  a  submarine  undergoing  similar 
maneuvers.  One  remaining  task  would  be  to  perform  maneuvering  computations  for 
which  experimental  data  exists  for  comparison. 


Figure  1.  SUBOFF  with  four  stern  appendages 
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Figure  2.  Velocity  profile  on  the  hull  for  CH-grid;  Ret  =  1,2  x  10 
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Figure  3.  Velocity  profile  on  the  hull  for  CH-grid;  Rol  -  1«2  x  10^ 


SUBOFF  WITH  STERN  APPENDAGES,  U+  vs  Y+  on  the  fins 


Figure  4.  Velocity  profile  on  a  stern  appendage  for  CH-grid;  Rcl  “  1.2  x  10 


SUBOFF  WITH  STERN  APPENDAGES  AT  0.0  DEGREE  ANGLE  OF  DRIFT 


Figure  5.  Convergence  history  for  axial  force  coefficient;  =  Rcl  1.2  x  10 
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Figure  6.  Convergence  history  for  axial  force  coefficient;  =  Ret  1.2  x  10 
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Figure  8.  Block  structure  in  the  stern  region: 


Figure  9.  Grid  around  a  stern  appendage: 
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Figure  11.  Velocity  profile  for  OH-grid;  Rol  =  1.2  x  10 
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igiire  12.  Convergence  history  comparison  for  axial  force  coefficient- 


O - O  Experimental  data 


efficient  versus  angle  of  d: 


SUBOFF  Bare  Flull  +  4  Stern  Appendages  at  O'"  Angle  of  Attack  (OH  grid) 
Pressure  Coefficient  on  the  Stern  Appendages 


Figure  14.  Pressure  coefficient  on  the  stern  appendages  for  the  OH-grid 


SUBOFF  Bare  Flull  +  4  Stern  Appendages  at  0®  Angle  of  Attack  (CFl-grid) 
Pressure  Coefficient  on  the  Stern  Appendages 


Figure  15,  Pressure  coefficient  for  the  CH-grid 


SUBOFF  Bare  Hull  +  4  Stern  Appendages  at  0°  Angle  of  Attack  (OH  grid) 
Skin  Friction  Coefficient  on  the  Stern  Appendages 
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Figure  16.  Sldn  friction  on  the  stern  appendages  for  OH-grid 


Time  history  of  the  rising  maneuver 

Z-force  on  top  surface  of  each  stem  plane  and  vehicle  inertial  position 
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e  distribution  during  the  pitch-up  maneuver 


Time  history  of  the  rising  maneuver 

Z-force  on  top  surface  of  each  stern  plane  and  vehicle  rotation  rates 


Figure  24 


Time  history  of  the  rising  maneuver 

Z-force  on  top  surface  of  each  stem  plane  and  vehicle  inertial  orientation 
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Time  history  of  the  diving  maneuver 

Z-force  on  bottom  surface  of  each  stem  plane  and  vehicle  inertial  position 


Figure  26. 


Time  history  of  the  diving  maneuver 

Z-force  on  bottom  surface  of  each  stem  plane  and  vehicle  rotation  rates 
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Figure  27. 


Time  history  of  the  diving  maneuver 

Z-force  on  bottom  surface  of  each  stem  plane  and  vehicle  inertial  orientation 
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APPENDIX  A 


NLIMERICAL  PREDICTIONS  OETURBLUENT  FLOWS  OVER 
UNDERWATER  VEHICLES  AT  HIGH  INCIDENCE 


Farhad  Davoudzadeh,  David  Boger,  Howard  Gibeling 
Computational  Mechanics  Department 
Applied  Research  Laboratory 
The  Pennsylvania  State  University 
State  College,  PA  16804 


The  high  Re>iiolds  number  turbulent  flow  over  an  appended  submarine 
geometry  at  high  incidence  angles  is  predicted  by  solving  the  governing 
Revnolds-Averaged  Navier-Stokes  equations.  These  flowfields  are  dominated 
by  leeward-side  vortical  structures  and  regions  of  crossflow  separation.  An 
algebraic  turbulence  model  as  w^ell  as  various  modified  versions  and  a  two- 
equation  model  are  used  in  several  sets  of  computations.  The  algebraic  model 
was  modified  to  account  for  the  effects  of  the  leeward-side  vortices  on  the 
underlying  viscous  layers.  The  computed  results  were  compared  with 
experimental  data,  and  advantages  and  disadvantages  of  the  models  are 
highlighted. 


INTRODUCTION 

Underwater  vehicles  and  weapons  are  required  to  be  highly  maneuverable.  This  means  flight  at  high 
incidence  angles.  Furthermo'e,  a  detailed  description  of  the  high  incidence  flowfield  is  needed  by  the 
designers  of  the  body  hull  hull  appendages  and  propulsors.  In  addition,  numerical  prediction  of  the 
trajectory  of  the  maneuvering  bodies  requires  computations  of  tlie  high  incidence  flowfield.  McDonald  and 
Whitfield  [1]  summarized  the  development  of  a  new  physics-based  method  for  the  prediction  of  the 
trajectory  of  the  underwater  vehicles  by  coupling  of  tlie  fluid  dynamics  with  the  vehicle  dynamics.  The 
equations  governing  the  vehicle  motion  are  commonly  referred  to  as  the  6DOF  equatiais.  Tlie  6DOF 
equations  describe  the  acceleration  of  a  rigid  body  given  the  totality  of  forces  and  moments  acting  on  it. 
The  RANS  equations  provide  the  means  to  calculate  the  hydrod>'namic  component  of  these  forces  and 
moments.  Integration  of  the  6DOF  equations  yields  the  time  history  of  the  vehicle’s  velocity  and  rotation 
rate.  From  this  information,  the  vehicle’s  trajectory  (the  history  of  tlie  position  and  orientation)  can  be 
deduced  by  integration  of  purely  kinematic  relations.  Davoudzadeh  and  others  [2]  implemented  the  method 
and  computed  trajectories  with  the  coupled  solvers  for  a  submarine  geometry  with  a  five-bladed  propeller. 
It  is  clear  that  the  accurate  prediction  of  the  forces  and  moments  is  key  to  the  accurate  prediction  of  the 
vehicle  trajectory.  Tlierefore,  tliis  paper  focuses  on  the  accurate  prediction  of  forces  and  moments  acting 
on  underwater  vehicles  at  high  incidence,  and  all  of  the  comparisons  made  between  the  computed  and  the 
experimental  data  are  for  the  force  components. 

The  flovrTields  associated  with  bodies  at  high  incidence  are  three-dimensional,  turbulent  and  contain 
regions  of  crossflow  separation  accompanied  by  leeward-side  vortical  structures.  Numerical  predictions  of 
these  flowfields  require  a  very  fine  grid  to  capture  the  off-body  vortical  structures  and  to  resolve  the 
boundary  layer.  The  ofl'-body  vortical  structures  are  essentially  inviscid  and  are  aeated  from  the 
underlying  turbulent  viscous  boundary  layer.  Therefore,  the  choice  of  the  turbulence  model  plays  a  major 
role  on  the  development,  structure,  stretching  and  strength  of  the  vortices.  Degani  and  others  [3] 
introduced  a  modification  to  the  widely  used  Baldwin-Lomax  eddy-viscosity  turbulence  model  that 
accounts  for  the  presence  of  the  leeward-side  vortex  structure.  They  presented  results  for  supersonic  [4] 
and  subsonic  [5]  high-incidence  turbulent  flow  obtained  using  a  parabolized  Navier-Stokes  method  [6]  and 
their  modified  turbulence  model.  They  reported  good  comparisons  with  the  experimental  data. 


Ib  the  current  work  the  models  examined  are  the  original  Baldvvin-Lomax[7],  Degani’s  modified  version  of 
the  Baldwiii-Lomax[3],  Baldwin-Lomax  with  various  search  parameters,  and  a  new  version  of  the  model 
which  has  m(xe  generality.  We  also  include  results  using  a  two-equation  turbulence  model  for  comparison. 
All  of  the  models  are  incorpo'ated  in  the  incompressible  Reynolds-Averaged  Navier-Stokes  (RANS)  code 
described  beIow\ 

Implies^  Unsteady  Incompressible  Flow  Solver 

The  three-dimensional  unsteady  incompressible  Reynolds- Averaged  Navier-Stokes  equations  are  first 
written  for  a  time-dependent  curvilinear  coordinate  system.  An  artificial  time  term  based  on  the  artificial 
compressibility  formulation  of  Chorin  [8]  is  then  introduced  to  facilitate  subiteration  at  each  physical 
time  step  (e.g.,  [9],  and  [10]),  Tlie  governing  equations  are  given  by 

T  A 
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+  k^p 


d  ^  +  k^u  +  kyV  +  k.w' 

K  :::  =  U  fOtk  =  ^ 

::=.V  for  k  ^  q 
K=  =  W  for  C 

In  these  equations,  p  is  static  pressure;  u,  v,  and  w  are  the  velocity  components  in  Cartesian  coordinates 
y,  and  z;  U,  V,  and  W  are  the  contravariant  velocity  components  in  the  curvilinear  coordinate  directions 
q,  and  ^  respectively;  and  x  is  time.  The  terms  T^.  ,7^  with  k  =  q,  or  are  the  viscous  flux 
components  in  curvilinear  coordinates.  The  Jacobian  of  the  inverse  transformation  is  J,  and  k^,  ky,  k^,  and 
ki  with  k  =  q,  or  are  the  transformation  metric  quantities,  where  a  subscript  denotes  differentiation. 
The  artificial  compressibility  parameter  [3  typically  has  a  value  of  5  to  10.  Finally,  I-,:  =  diag  [0,1, 1,1],  and 
x’  is  an  artificial  time  associated  with  the  subiteration.  The  first  term  in  Eq.  (1)  vanishes  as  the  subiteration 
converges,  giving  a  solution  to  the  physical  unsteady  equations. 

A  thin-layer  approximation  is  used  to  simplify  the  viscous  stress  terms,  and  several  turbulence  models  have 
been  implemented  in  modular  form,  including  algebraic  turbulence  models  [7],  and  a  q-co  model  [11]. 
Details  of  treating  the  viscous  terms  are  given  by  Gatlin  [12],  with  an  improvement  due  to  Chen  [13] 
regarding  the  computation  of  tangential  velocity  derivatives  normal  to  a  solid  surface.  This  simple 
technique  works  extremely  well  on  highly  skewed  grids  [13].  - 


Equation  (1)  is  first  discretized  as  a  cell-centered  finite-volume  approximation,  and  the  cell-face  numerical 
flux  vectors  are  obtained  using  Roe's  [14]  flux  difference  split  scheme.  The  nonlimited  form  of  the 
dependent  variable  extrapolation  method  of  Anderson,  Thomas,  and  van  Leer  [15]  is  used  to  obtain  second 
and  third-order  flux  vectors  [15]-[16].  The  third-order  upwind-biased  numerical  flux  [15]  is  used  for  the 
present  results. 

The  implicit  solution  procedure  for  tlie  disaetized  equations  varies  according  to  whether  the  solution  of 
interest  is  steady  or  unsteady.  For  steady  flows,  the  second  term  in  Eq.  (1)  is  omitted,  reducing  the 
equations  to  a  standard  artificial  compressibility  formulation.  The  spatial  difference  terms  are  linearized 
to  produce  a  linear  system  for  A(y,  which  is  solved  iteratively  at  each  time  step  [17].  In  tins  case,  both 
AQ^  and  the  steady  spatial  derivative  residuals  are  zero  upon  convergence  of  the  (outer)  time-step  iteration 
to  a  steady  solution.  For  unsteady  flows,  an  unsteady  residual  equation  is  defined  from  the  last  four  terms 
in  Eq.  (1)  and  is  solved  iteratively  for  In  one  dimension,  for  example,  the  unsteady  residual  becomes 


Fcg"^')  =  /  =  0 

At 


The  function  9(x)  in  (2)  is  solved  by  an  approximate  Newton  method  given  by 


(2) 


n+l,/n 


(3) 

where  m  =  1,  2,  3, ...  is  the  Newton  iteration  index,  and  9  '(x)  is  the  Jacobian  matrix  of  the  vector 
9(x)  modified  by  replacing  T  by  I  to  utilize  artificial  compressibility.  In  principle,  the  generated  sequence 
Qo+i,  m+i  2a\d  hence  Eq.  (2)  is  satisfied.  Tlie  linear  system  for  in  Eq.  (3)  is 

solved  iteratively  for  each  Newton  iteration. 


A  s>mmeLric  Gauss-Seidel  relaxation  is  used  to  solve  tiie  linear  system  of  equations  associated  with  each 
iteration  of  Eq.  (3).  Because  the  flux  Jacobian  of  the  flux  vector  based  on  Roe’s  formulation  is  difficult  to 
obtain  anal>tically  in  three  dimensions,  and  also  in  the  interest  of  simplicity,  the  flux  Jacobian  is  computed 
numerically  [18]  [19].  This  solution  scheme  is  referred  to  as  a  discretized  Newton-relaxation  [20],  or  the 
DNR  scheme  [19].  A  multigrid  scheme  [21]  is  used  to  accelerate  the  numeriail  solutions,  and  this  scheme 
has  been  extended  to  multiblock  [22]  and  unsteady  flows  [23].  The  present  solutioi  procedure  is,  therefore, 
a  muitigrid  scheme  for  three-dimensional  unsteady  viscous  flow.  This  procedure  is  the  basis  for  a  flow 
simulation  code  called  UNCLE  (UNsteady  Computation  of  fieLd  Equations). 


Turbulence  Models 


AJgebraic  Models 

Tlie  well-known  Baldwin-Lomax  turbulence  model,  which  is  based  on  tlie  two-layer  model  of  Cebeci  and 
others  [24],  does  not  take  into  account  the  complexities  of  the  tliree-dimensional  flowfield.  The  model  gets 
particularly  confused  when  local  vorticity  exists  in  areas  other  than  the  near-wall  boundary  layer.  This 
situation  happens  when  the  flowfield  contains  vortical  structures  off  of  the  body.  This  is  due  to  the  fact  that 
tlie  eddy  viscosity  in  the  outer  region  is  a  function  of  the  magnitude  of  tlie  local  vorticity.  In  the  original 
Baldwin-Lomax  model  formulation  the  algebraic  eddy  viscosity  jit  is  given  by: 

n  =l(M,),nner>  ^  ^  Vc 

\(M,)cu,er^  y>yc 


(4) 


where  y  is  the  normal  distance  from  the  wail  and  yc  is  the  smallest  value  of  y  at  which  values  from  the  inner 
and  outer  formulas  are  equal.  The  Prandtl-Van  Driest  formulation  is  used  in  the  inner  region 


where 


l(ol  is  the  magnitude  of  the  vorticity 


(5) 


(6) 


co:  +  Cl);  +  co; 


(7) 


and 


=i.4p.^JpJy 


(8) 


[n  the  outer  region,  for  attached  boundary  layers,  the  turbulent  viscosity  coefficient  is  given  by 
(MXu,er  =  ^^cpP  F.-o>ceFHeb  (v) 


(9) 


where  K  is  the  Clauser  constant,  Cep  is  an  additional  constant  and  Ftteb  is  the  Klebanoff  intermittency 
factor,  and 


F 


wake 


F 


max  ■ 


The  quantities  Y„iax  and  Fmax  are  derived  from  the  function 


(10) 


F(y)  =  |ty| 

(11) 

where  Fma.x  is  the  maximum  value  that  the  function  F(y)  takes  in  a  local  profile,  and  Ymax  is  the  value  at 
which  Fma.x  occurs.  The  constants  appearing  in  the  above  relationships  have  been  determined  by  requiring 
agreement  with  the  Cebeci  [24]  formulation  for  constant  pressure  boundary  layers  at  transonic  speeds.  The 
values  determined  are  K  =  0,0168,  A'^  =  26,  and  Cep  =1.6. 

As  previously  demonstrated  by  Degani  and  Schiff  [3]  the  main  difficulty  in  applying  the  Baldwin-Lomax 
turbulence  model  to  three-dimensic^al  flow  containing  regions  of  crossflov/  separation  is  the  determination 
of  the  length  scale  Ymax  and  in  turn,  of  evaluating  (iiJouter  for  boundary-layer  profiles  in  the  crossflow 
separatiCHi  region.  To  demcxistrate  this  difficulty,  they  looked  at  behavior  of  the  F(y)  function  along  two 
ra>'s-one  located  in  the  windward  side  at  (})  =  <{)t  and  the  other  one  on  the  leeward  side  at  <{)  =  shown  in 
Fig.  1.  The  functions  are  shown  schematically  in  Figs.  2a  and  2b,  respectively.  The  function  along  the  ray 
in  the  windward  side  contains  only  one  well-defined  peak,  whereas  along  the  ray  in  the  leeward  side  it 
contains  two  distinct  peaks.  These  peaks  are  manifestations  of  the  vorticity  magnitude  contained  in  the 
function  F  (y).  The  first  peak  in  Fig.  2b  comes  from  the  vorticity  magnitude  contained  in  the  boundary 


layer;  whereas  the  second  peak  onginates  from  the  vorticity  magnitude  contained  in  the  off-wal’l  vortical 
structure.  If  one  uses  the  length  scale  corresponding  to  this  second  peak,  as  would  be  tlie  case  for  the 
original  Baldwin-Lomajc  model,  the  computed  eddy  viscosity  in  Eq.  (9)  will  be  at  least  an  order  of 
magnitude  larger  than  the  eddy  viscosity  value  related  to  the  length  scale  corresponding  to  the  first  peak. 
This  exaggerated  eddy  viscosity  will  have  the  eftect  of  weakening,  deforming,  and  dissipating,  if  not 
eliminating  all  together,  the  vortical  structures  in  the  flowfield.  Furthermore,  the  computed  drag  on  the 
body  will  be  too  high.  This  has  serious  implications  on  the  prediction  of  the  maneuverability  of  the 
underwater  vehicles  with  such  a  model. 


Fig  1.  Schematic  of  Flow  Structures  in  Crossflow  Plane  (After  Degani  etrd.  [3]) 
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Fig. 2  Behavior  of  F(y)  at  Large  Incidence  at  a  Cross  Section  Containing  Vortices  (After  Degani  el  al.  [3]) 


Partial  Differential  Turbulence  Morleis 

The  idea  Ixthind  eddy  viscosity  mtxlels  comes  from  the  Boussinesq  approximation  that  the  Reynolds  shear 
stress  can  be  related  to  the  mean  shear  by 

-liv  =  V^dU  /3v. 

(12) 

Equation  (12)  defines  tlie  turbulent  viscosity,  Vj,  which  is  dimensionally  equal  to  die  product  of  a  velocit}^ 
scale  and  an  appropriate  lengdi  scale.  For  example,  die  Baldwin-Lomax  model  described  above  chooses 
the  velocity  scale  to  be  the  magnitude  of  vorticity  times  a  length  scale,  and  the  length  scale  is  prescribed 
algebraically.  Unfortunately,  prescribing  a  length  scale  while  maintaining  sufficient  generality  is  difficult, 
as  was  just  described. 


A  method  which  provides  for  more  generality  is  to  incorporate  partial  differential  equations  for  the 
transport,  production,  and  dissipation  of  botli  the  velocity-  and  iengtii-scale  determining  quantities.  Here, 
following  Coakley  [11],  the  velocity  scale  is  taken  to  be  the  square  root  of  turbulent  kinetic  energy,  q,  and 
the  length-scale  determining  quantity  is  the  time  rate  of  dissipation,  co,  governed  by 
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The  eddy  viscosity  is  then  formulated  as 


(13) 


Mt  =  pCuDq'  /CO 


(14) 


where  D=i-exp(-a4v/v)  is  a  damping  function  based  on  y,  tlie  distance  to  the  nearest  solid  wail. 
Completing  tlie  model,  a=0.020,  Ci=0.055-f-0.5D,  C2=  0.833,  Cu  =  0,09.  raui  Oq  =  =  2.0. 

The  solution  metliod  is  tlie  same  as  tlie  method  described  previously  for  tlie  flow  equations,  except  tliat  it  is 
applied  to  two  equations  instead  of  four.  The  velocity  field  is  treated  as  known  for  tliis  part  of  the  solution, 
but  the  turbulence  equations  are  coupled  together  tlirough  tiieii*  source  terms.  More  details  of  the  mo<jei 
and  its  implementation  are  described  elsewhere  [25]. 


RESULTS 


A[[  of  tlie  computations  presented  here  are  conducted  on  a  SUBOFF  body  with  four  stern  appendages. 
SUBOFF  has  tlie  shape  of  a  typical  underwater  vehicle.  Groves  and  otliers  [26]  designed  an  axisymnietric 
hull  model  for  SUBOFF  that  could  be  used  with  mid  without  ivpical  appendage  components,  such  as  a 
fairwater  (sail)  and  four  identical  stem  appendages.  Rtxldy  [27]  reported  experiments  tliat  include 
measuremenls  of  forces  and  moments  acting  on  die  SUBOFF  body,  with  mid  without  appendages  at  drill 
angles  ranging  from  0.0°  to  18.0°.  All  experiments  were  airried  out  at  a  Reynolds  number  of  14.000,000 
based  on  tlie  bo/dy  iengtli.  Since  the  body  is  symmetric,  only  half  of  tlie  body  is  gridded  for  die 
computations.  A  C-type  grid  widi  8  blocks  and  81x73x33  (195, 129)  grid  points  per  block-giving  a  total  of 
1,561,032  grid  points-  is  used.  The  grid  is  split  into  two  pjirts  in  die  axial  direction  mid  four  parts  in  die 
circumferential  direction.  No  blocking  is  made  in  die  radial  direction.  Each  fin  surface  has  17  grid  points 
in  die  axial  direction  and  33  grid  points  in  die  radial  direction. 

Original  Baldwin-Lomax 

Using  the  algebraic  model  as  proposed  by  Baldwiii-Loraax  [7],  a  series  of  steady  computations  for  different 
drift  migles  are  performed.  Fig.  3  shows  contours  of  axial  components  of  vordcity  on  SUBOFF  at  18° 
angle-of-drift  for  the  computadons  using  the  original  Baldwin-Lomax  model.  Tliis  dgure  shows  die  diree 
dimensionality  of  the  flow,  and  the  development  and  roll-up  of  the  vortex  structure  on  the  leeward  side  of 
the  body.  However,  as  pointed  out  earher,  the  high  value  of  the  calculated  eddy  '/iscosity  associated  with 
the  second  peak  (see  Fig.  2b)  causes  the  vortex  to  be  w^eak  and  not  well-defined.  Fig.  4  shows  a  cross 
section  of  die  computed  pressure  distribution  at  =  0.63.  Fig.  5  shows  the  axial  force  against  the  drift 
migle,  witli  die  error  bars  associated  with  the  experimental  data.  The  computed  axial  force,  for  most  angles 
of  drift,  underestimates  the  experimental  values,  i.e.,  extra  eddy  viscosity  induces  higher  drag  ou  die  body. 
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Fig.  3  Contours  of  Axial  Component  of  Vorticity 
on  SUBOFF  (Original  Baldwin-Lomax,  (3=18  deg, 
Re  =  14M) 
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Fig.  4.  Contours  of  Pressure  in  the  Cross 
Section  at  xyL=0.63  (drift  angle  =18  deg., 
Re=14M;  Original  Baldv/in_Lomax 


Fig.  5  Computed  and  Measured  Axial  Force  Coefficient  for  Various  ^Angles  of  Drift  (Original  Baldwin- 
Lomax) 


Baldvvin-Lomax  with  constant  search  parameter 

In  die  second  series  of  computations,  tlie  search  for  along  radial  rays  (shown  in  Fig.4)  is  limited  to  liie 
value  y'^  =  750.  Once  tliis  Fnnu  is  found,  tlie  corresponding  Ymax  is  detennined.  Eddy  viscosity  is 
subsequently  determined  from  Eqs.  9  tmd  10.  Figure  6  shows  tiie  computed  contours  of  the  x-compoaent 
of  vorticity  for  18.0°  angle-of-drift  computations.  Comparison  of  tins  figure  widi  Fig.  3  demonstrates  die 
formation  of  a  stronger  and  better  defined  leeward  vortical  structure.  A  cross  section  of  pressure 
distribution  at  x/1  =  0.63  is  shown  in  Fig.  7.  hi  contrast  widi  the  Fig.  4,  the  presence  ot'  a  pair  of  very  well 
defined  vortices  is  evident  in  the  leeward  side  of  die  body.  Figure  7  also  shows  die  radicil  dii*ections  along 
which  the  function  F  is  evaluated.  The  radial  lines  of  161.26°  and  167.08°  pass  almost  dirough  die 
pmnary  vortex  core.  .A  secondary  vortex  also  occurs  between  die  radial  lines  of  128.8°  and  145.76°. 
Variation  of  F  against  yL  along  the  radial  lines  shown  in  Fig.  7,  for  18°  drift  angle  is  presented  in  Fig.  8b. 

In  diis  figure,  die  values  of  F,uax  picked  by  die  computer  code,  based  on  die  search  pai'ameter  of  750  are 
also  shown.  Figure  8b  shows  that  for  the  most  part,  F^u^x  is  related  to  die  boundary  layer  diickiiess  -  i.e.  the 
first  peak  is  realized  -  except  at  145.76°.  The  lower  value  of  F^ax  calculated  in  diis  computation  results  in 
lower  values  for  the  eddy  viscosity  relative  to  the  original  Baldwin-Lommx  model  As  a  result  the  leewai'd 
vortex  is  better  defined  and  maintained.  Figure  9  shows  comparisons  of  the  computed  and  experimental 
data  for  various  angles  of  drift.  The  model  prediction  compares  well  widi  all  the  positive  angles  of  drift 
higher  than  10°.  Tlie  model  overpredicts  the  axial  force  for  drift  angles,  between  0°  and  10°.  An 
examination  of  the  function  F  along  the  zero  radial  line  for  the  zero  angle-of-attack  cases  shown  in  Fig.  8a, 
reveals  that  due  to  the  search  constraint  of  750,  a  lower  value  of  F^ax  is  picked  by  the  code.  This  explains 
why  the  magnitude  of  axial  force  in  Fig.  9,  for  small  drift  angles  is  overpredicted.  Obviously  for  small  drift 
angles,  one  could  choose  a  higher  value  of  the  search  parameter,  but  for  consistency,  ail  computations 
presented  in  Fig.  9  are  conducted  with  the  search  value  of  750. 
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Fig.  6  Contours  of  Axial  Component  Figure  7.  Contours  of  E^essure  in  the  Cross  Section  at 

of  Vorticity  on  SUBOFF  (Search  Parameter  x/L=0.63  (drift  angle  =  18  deg.  Re  =  14M; 

=  750,  P  =  1  8  deg.  Re  =  14M)  Baidwin-Lomax  witii  Constant  Search  Parameter  of 


750) 


Fig.  8  The  F  as  Functions  of  Above  die  Body  Fig  9  Computed  mid  Measured  Axial  Force 

Surface  (Semxh  Parcuneter  =  750)  Coefficient  for  Vmious  Angies  of  Drift 

(Search  Parameter  =  750) 

A  close  ex^miination  of  Figs.  8a”8b  reveals  Quit  a  search  pm'ameter  of  -  2,500  is  most  likely  to  capture 
Qie  first  peak  for  all  drift  angles  under  consideration.  Fig.  10  shows  Qie  comparison  between  Qie  computed 
axial  force  coefficient  obtained  wiQi  search  ptirameter  of  2,500  and  Qie  experimental  values  at  various 
angles  of  drift.  For  drift  angles  higher  than  approximately  10°,  Qie  comparison  is  similar  to  Qie  results 
obtained  wiQi  the  search  parameter  of  750.  However,  for  all  Qie  angles  between  0°  and  10°,  the  prediction 
is  remarkable.  Unfortunately,  because  of  tlie  asymmetry  of  tlie  experimental  data  (for  a  symmetric  body), 
no  model  can  predict  these  experimental  data  for  both  negative  and  positive  values  of  Qie  drift  angle. 
Comparison  of  the  lateral  force  coefficient  is  shown  in  Fig.  11.  It  should  be  noted  Qiat  the  lateral  force 
coefficients  compare  well  with  the  experimental  data  for  all  of  Qie  models  discussed  above.  This  is  due  to 
the  fact  that  the  contribution  of  the  viscous  stresses  to  the  computed  component  of  the  force  is  much  higher 
in  the  axial  direction  Qian  in  the  lateral  direction.  The  lateral  force  comes  mainly  from  the  pressure  field. 
This  is  better  demonstrated  in  Figs,  12  and  13.  Figure  12  shows  the  contributions  of  the  viscous  stresses 
and  pressure  field  to  the  total  axial  force  coefficient.  Figure  12  also  shows  the  total  axial  force  coefficient, 
both  experimental  and  computed.  For  drift  angles  greater  than  10°,  the  pressure  field  contribution  to  Qie 
axial  force  creates  a  Qinist  on  the  body  (due  to  the  shape  anomaly),  whereas  for  drift  angles  less  than  10°, 
the  contribution  is  drag  on  the  body.  Contribution  of  the  viscous  stresses  to  the  axial  force  is  ah  drag.  Fig. 
13  shows  similar  divisions  for  the  lateral  force  coefficient.  It  is  evident  from  Qiis  figure  that  the 
contribution  of  the  viscous  stresses  to  the  lateral  force  is  very  small,  compared  to  pressure. 


Axial  force  coefficient  Axial  force  coefficient 


Fig.  10  Computed  and  Measured  Axial  Force 
Coefficient  for  Various  Angles  of  Drift  (Search 
Parameter  -  2,500) 


Angle  of  drift 


Fig.  11  Computed  &  Measured  Lateral  Force 
Coefficieiu  for  Various  Mgles  of  Drift 
(Search  Parameter  =  2,500) 


Fig.  12  Computed  Contribution  of  Viscous  Stresses  Fig.  13  Computed  Contribution  of  Viscous 

&  Pressure  to  tlie  Axial  Force  Coefficient  Stresses  and  Pressure  Field  to  die  Lateral 

Force  Coefficient 


Towards  a  generalized  model 

The  modified  Baldwia-Lom.ax  witli  die  constant  search  parameter  of  =  2,500  produces  results  drat 
compare  very  well  wtdi  die  experimental  data.  Choosing  an  appropriate  search  parameter  in  this  way  relies 
on  trial  and  error  and  on  runnmg  one  or  two  cases,  preferably  at  O.O'^  and  a  high  degree  incidence  angle. 
Degani  and  Schiff  [3]  suggested  other  strategies  for  setting  a  search  limit  which  require  die  prior 
loiowledge  of  the  windward  side  and  the  leew'ard  side  of  the  body.  Tlie  issue  widi  maneuvering  bodies  is 
that  the  windward  and  the  leeward  sides  are  not  fixed.  The  vehicle,  in  the  course  of  a  maneuver,  causing  a 
significant  change  its  orientation,  causing  a  change  in  the  kxations  of  the  leeward  and  windward  sides. 

To  eliminate  the  need  for  a  constant  search  parameter  and  the  need  to  know  the  leeward  and  the  windward 
sides  of  the  vehicle  a  prion,  a  new  model  is  suggested.  The  e-ssence  of  this  new  model  is  that,  in  the 
boundary  layer,  tiie  magnitude  of  die  axial  component  of  vorticity  is  very  small,  compared  widi  diat  of  the 
lateral  component  of  vorticity.  Therefore,  ©x  tloes  not  contribute  much  to  the  magnitude  of  the  vorticity  in 
the  boundary  layer.  However,  outside  of  the  boundary  layer,  where  the  off-body  vortex  structures  exist,  the 


contribistion  of  the  co^  to  the  voricity  magnitude  is  large  (see  Fig.  2b).  This  implies  that  cOx  can  be  used  to 
eliminate  the  second  peak  in  F  by  defining  a  new  function  G  as  follows: 


where, 


(15) 
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(16) 


The  functions  F  and  G  are  shown  schematically  in  Fig.  14.  The  only  difference  between  the  two  functions 
is  the  way  in  which  I  ST  I  is  defined.  It  is  assumed,  as  shown  in  Fig.  i4a,  that  G  eliminates  the  second, 


larger  peak  produced  by  F,  and  that  the  peak  in  G  occurs  at  approximately  the  same  y  location  as  the  first 


peak  in  F.  This  is  a  fairly  sound  assumption,  since  the  contribution  of  cOx  in  the  boundary  layer  to  the 
vorticity  magnitude,  and  tlierefore  die  size  of  die  peal:,  is  negligible.  Furthermore,  since  the  contribution 
of  oOx,  outside  of  the  boundary  layer  where  vortices  exist  is  large,  its  exclusion  from  Eq.  (16)  can  subdue  the 
large  peak  caused  by  the  vortical  structures.  For  IS.O""  angle-ol- incidence  compulations,  Fig.  15  shows  the 
variation  of  different  forms  of  F  and  G  against  die  radial  distance  from  die  wail  at  (j)  =  IbT.OS'^.  Figures 
15a- 15b  and  15c  shows  F  when  only  the  cOx,  cOy  or  cOz  were  considered  in  Eq,  (7).  Figures  I5d  md  15e 
show  functions  G  and  F.  It  is  clear  from  Fig.  15a  that  the  magnitude  of  COx  is  not  noticeable  in  die  boundai'y 
layer  region,  and  it  is  very  large  in  the  vortex  region.  Comparison  of  Fig.  I5a  with  Fig.  15e  demonstrates 
how  diis  large  value  of  cOx  has  dominated  the  overall  magnitude  of  die  vorticity  and  has  produced  the  Icirgei' 
peak  in  F.  On  die  odier  hand.  Fig.  15d  shows  how  die  search  tbr  die  maximum  peak  in  G  fields  a  length 
scale  dial  is  (he  same  as  die  first  peak  in  F,  Fig.  15e.  Using  diis  new  modified  mediod  a  series  of 
computations  is  performed  for  incidence  angles  of  0.0°  through  18.0°.  Figure  16  shows  the  pressure 
distribution  at  x/1  =  0.63  for  the  drift  angle  of  18.0°.  The  leeward  vortices  are  clccirly  shown  on  the  top 
portion  of  die  figure.  These  vortices  are  stronger  and  better  defined  dum  diose  obtained  under  die  original 
Baldwiii-Lomax  model,  but  they  are  not  as  strong  as  die  ones  computed  with  die  constant  search 
parameters.  An  examination  of  the  variation  of  F  along  radial  lines,  Fig.  17,  shows  diat  aldiougii  die  model 
picks  die  appropriate  peaks  for  die  determination  of  Fhui.k  for  most  rays,  it  fails  for  some,  e.g.  for  ^  ~ 
174,87°.  This  deviation  occurs  because  the  vortical  structures  ai'e  not  aligned  with  the  axial  direction  of  die 
body,  causing  die  calculated  lateral  vorticity  to  have  a  larger  value  tiian  die  value  if  die  vortex  core  was 
aligned  v\idi  die  axis  of  the  body.  This  causes  die  shape  of  G  to  be  different  than  the  assumptions 
presented  earlier.  Comparison  of  die  axial  force  coeffrcieiit  widi  die  experimental  data  for  diis  new  model 
is  shown  in  Fig.  18.  Unlike  die  odier  computalions,  die  computed  values  fail  widiin  die  error  bars  for  all 
positive  and  negative  angles  of  drift. 
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Figs.  14  Behavior  of  F(x)  &  G(y)  . 


Figs.  15a  through  15e  Variation  of 
Different  Forms  of  F  and  G  Against 
the  Radial  Distance  from  the  Wail 
at  0=167.08° 
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Fig.  16  Contours  of  E^essure  in  the  Cross  Section 
at  x/L=0.63  (Drift  Angie=18^  Re=14M,  New 
Model) 
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Fig.  18  Computed  and  Measured  Axial  Force  Fig.  17  The  F  as  Functions  of  Radial  Distance  Above  the  Body 

Coefficient  for  Various  Angles  of  Drift  (New  Surface  (New  Model) 

Model; 


Model 

Algebraic  models  are  highly  desirable  as  turbulence  models  because  they  are  inexpensive  computationally 
and  are  generally  very  robust.  On  the  other  hand  higher-order  models  may  be  required  to  attain  sufficient 
generality  and  accuracy.  Here,  Coakley’s  q-co  model  vv^s  applied  to  tlie  Fnest  grid  for  the  cases  of  0*^,  8°, 
and  18°  angle-oFattack.  The  lateral  force  was  well  predicted  similar  to  the  one  shown  in  Figure  11; 
however,  the  axial  force  compared  in  Fig.  19,  is  showm  to  be-  too  high.  It  was  suspected  that  the  oiuse  was 
an  over-prediction  of  turbulent  kinetic  energy  at  stagnation  points  (principally  at  the  nose).  A  treatment  for 
this  problem  has  been  suggested  by  Durbin  [28],  it  amounts  to  the  following  specified  on  co: 

co>^3i2S,jS,j)C^D 

(17) 

The  results  at  0°and  8°  are  somewhat  improved  by  this  modification,  see  Fig,  19,  but  at  the  cost  of 
robustness  -  -  the  case  at  18°  angle-of-attack  could  not  be  completed. 

For  both  twivequation  models,  the  cx»mputed  eddy  viscosity  is  relatively  high  in  the  vortex  cores,  and 
consequently,  the  shed  vorticity  is  diffused  too  quickly. 


Fig.  19  Computed  and  Measured  .4xial  Force  Coefficient  for 
Various  Angles  of  Drift  (Coakley’s  and  Diu-bin*s  q-co  Model) 

GRID  STUDY 

Previously  Ref.  [25]  investigated  the  influence  of  grid  density  on  tlie  accuracy  of  tlie  computed  forces  and 
moments,  witli  three  grid  were  constructed:  A  fme  grid  witli  12  blocks  and  81x89x41  grid  points  per  block, 
a  medium  grid  witli  8  b[cx:ks,  as  described  earlier,  and  a  coarse  grid  witli  8  blocks  and  41x37x17  grid  points 
per  block.  The  results  obtained  from  these  tliree  grids  considering  only  the  variation  of  the  axial  force 
coeffiicient  witli  angle  of  drift  showed  that  only  slight  differences  exist  between  die  fine  and  medium  grids. 
Computations  using  the  coarse  grid  under-predict  the  axial  force  coefficient  for  positive  angles  of  drift.  .All 
these  grids  compared  well  for  the  lateral  force  coefficient  and  yawing  moment  coefficient.  This  is  due  to 
die  fact  diat  the  contribution  of  the  viscous  stresses  to  the  lateralTbrce  and  yawing-moment  is  not  as  large 
as  it  is  for  die  axial  force. 


CONCLUSION 


A  series  of  computations  cire  performed  to  evaluate  die  suitability  of  several  idgebraic  as  well  as  partiai- 
diflerential  two-equation  turbulence  models.  The  algebraic  models  considered  ;;ire  the  original  Baldvvin- 
Lomcix  model,  Degani  mid  Schiff  variation  of  die  Bcildwin-Lomax  model,  and  a  new  propOvSed  mcxiel. 
Furdiermore,  the  two-equadon  model  of  Co^ikley  is  also  evaluated.  Tlie  computadons  demonstrate  that  die 
original  Baldwin -Lomax  is  not  suitable  for  die  diree-dimensional  computations  diat  include  vortical 
structure  away  from  die  Ixxly.  This  model  dissipates  and  distorts  die  vortical  structures  and  underesdmates 
axial  and  lateral  forces.  The  modified  Baiwin-Lomax  widi  die  specided  search  piirameters  prexiuces  very 
well-defined  vordces,  and  die  computed  forces  compare  well  widi  die  experimental  data  only  if  a  suitable 
search  parameter  is  used.  However,  die  search  for  the  suitable  search  pai'ameter  relies  on  (rial  and  error, 
and  dierefore,  diese  models  are  not  general.  The  proposed  algebraic  model  is  a  general  model  that  does  not 
need  a  search  parameter  and  does  not  require  prior  knowledge  of  die  leeward  or  the  windward  side.  It 
preserves  die  structure  of  die  vortices  belter  dian  die  original  Baidwin-Lomax  does.  However,  in  its  present 
form,  it  assumes  diat  die  core  of  die  vortical  structures  is  aligned  with  the  axial  axis  of  die  body.  A  more 
general  extension  of  die  model  should  consider  aligimient  of  the  vordcal  structures  in  any  direction.  The 
two-  equation  model  of  Coakley  is  also  considered.  The  model  predicts  relatively  high  eddy  viscosity  in  the 
cores  of  the  shed  vortices,  and  so  the  vortices  diffuse  too  quickly. 
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